Introduction to CoRC
Jonas Förster & Jürgen Pahle
2017-09-06
Case-studies and Workflows
Initial Setup
library(tidyverse)
library(parallel)
library(CoRC)
# helper to run tasks in parallel on all cores
mapInParallel <- function(data, fun, ..., .prep = {}) {
cl <- makeCluster(detectCores())
clusterCall(cl = cl, fun = eval, .prep, env = .GlobalEnv)
ret <- clusterApplyLB(cl = cl, x = parallel:::splitList(data, length(data)), fun = lapply, as_mapper(fun), ...)
stopCluster(cl)
flatten(ret)
}3D Trajectory Plot of a Calcium Model
This example loads the (Kummer2000 - Oscillations in Calcium Signalling) model. The model has 3 species which oscillate. These oscialltions can be visualized as a trajectory through a 3D space. The example does this once in a deterministic and once in a stochatic fashion.
loadSBML("http://www.ebi.ac.uk/biomodels-main/download?mid=BIOMD0000000329")
# Run 24 sec (2 Periods)
setTimeCourseSettings(24, intervals = 10000)
timeseries <- list(
deterministic = runTimeCourse()$result,
stochastic = runTimeCourse(method = list(method = "directMethod", use.random.seed = T, random.seed = 1))$result
)
# simplify species names so they are valid R syntax
timeseries <- map(timeseries, set_tidy_names, TRUE)
# Create the same plot for both timeseries
plots <- map(
timeseries,
plotly::plot_ly,
type = "scatter3d",
mode = "lines",
x = ~ G.alpha,
y = ~ activePLC,
z = ~ Calcium,
color = ~ Time
)
unloadModel()
plots$deterministicplots$stochasticStatistics of Repeated Stochastic Simulations
This implements an example from the (Condor-COPASI paper). The example illustrates advantages of parallel processing.
# Run 1000 stochastic time series possibly in parallel
loadModel("https://raw.githubusercontent.com/copasi/condor-copasi/master/examples/stochastic_test.cps")
# timeseries <- 1:1000 %>% map(~ runTimeCourse()$result)
timeseries <-
# Defines parallel evaluation:
mapInParallel(
# perpare worker (.prep),
.prep = quote({
library(CoRC)
loadModel("https://raw.githubusercontent.com/copasi/condor-copasi/master/examples/stochastic_test.cps")
setTimeCourseSettings(method = list(method = "directMethod", use.random.seed = T))
}),
# iteration data (1000 random seeds),
1:1000,
# iteration code (formula ~)
~ runTimeCourse(method = list(random.seed = .x))$result
)
# Combine all results and reshape the data
plotdata <-
timeseries %>%
bind_rows() %>%
group_by(Time) %>%
# calculate mean and sd for all time points
summarise_all(funs(mean, sd)) %>%
# gather all values so the column "name" identifies "a_mean", "b_sd" etc.
gather("name", "value", -Time) %>%
# split up information on species (a,b,c) and type of value (mean, sd)
separate(name, c("species", "type"), "_") %>%
spread(type, value)
print(plotdata, n = 6)
#> # A tibble: 2,403 x 4
#> Time species mean sd
#> * <dbl> <chr> <dbl> <dbl>
#> 1 0.00 a 8.003797 0.0000000
#> 2 0.00 b 8.003797 0.0000000
#> 3 0.00 c 8.003797 0.0000000
#> 4 0.05 a 7.064646 0.2432819
#> 5 0.05 b 8.114339 0.1177902
#> 6 0.05 c 5.596016 0.4555996
#> # ... with 2,397 more rows
plot <-
ggplot(data = plotdata, aes(x = Time, y = mean, group = species, tt_sd = sd)) +
geom_ribbon(aes(ymin = mean - sd, ymax = mean + sd, fill = species), alpha = 1 / 4) +
geom_line(aes(color = species)) +
guides(fill = "none") +
labs(
x = paste0("Time (", getTimeUnit(), ")"),
y = paste0("Concentration (", getQuantityUnit(), ")"),
color = "Species"
)
unloadModel()
plotly::ggplotly(plot, tooltip = c("group", "x", "y", "tt_sd"))Parameter Scan
2D Scan Over the Cartesian Product of Two Species Concentration Vectors
This implements an example from the (Mendes2009 paper) on COPASI use cases.
loadSBML("https://www.ebi.ac.uk/biomodels-main/download?mid=BIOMD0000000068")
setSteadyStateSettings(calculateJacobian = FALSE, performStabilityAnalysis = FALSE)
# Cartesian product of the input values
scan <- cross_df(
list(
cysteine = 0.3 * 10 ^ seq(0, 3, length.out = 6),
adomed = seq(0, 100, length.out = 51)
)
)
scan <-
scan %>%
mutate(
# Calculate steady state fluxes for every row
ss_fluxes = pmap(., function(cysteine, adomed) {
setSpecies("Cysteine", initial.concentration = cysteine)
setSpecies("adenosyl", initial.concentration = adomed)
ss <- runSteadyState()
stopifnot(ss$result == "found")
ss$reactions$concentration.flux
}),
# The second flux is CGS
CGS = map_dbl(ss_fluxes, 2),
# The third flux is TS
TS = map_dbl(ss_fluxes, 3)
)
plot <-
ggplot(data = scan, aes(x = adomed, group = cysteine)) +
geom_point(aes(y = CGS, color = "CGS")) +
geom_point(aes(y = TS, color = "TS")) +
geom_line(aes(y = CGS, color = "CGS")) +
geom_line(aes(y = TS, color = "TS")) +
labs(
x = paste0("Adomed (", getQuantityUnit(), ")"),
y = paste0("Flux (", getQuantityUnit(), " / ", getTimeUnit(), ")"),
color = "Species"
)
unloadModel()
plotly::ggplotly(plot)2D Scan over Random Concentrations of Two Species
This implements an example from the (Mendes2009 paper) on COPASI use cases. It is in many ways similar to the previous example but is written to run parallelized.
loadSBML("https://www.ebi.ac.uk/biomodels-main/download?mid=BIOMD0000000068")
# 10000 repeats of steady state task with random cysteine and adomed
scan <-
# Defines parallel evaluation:
mapInParallel(
# perpare worker (.prep),
.prep = quote({
library(CoRC)
loadSBML("https://www.ebi.ac.uk/biomodels-main/download?mid=BIOMD0000000068")
setSteadyStateSettings(calculateJacobian = FALSE, performStabilityAnalysis = FALSE)
}),
# iteration data (10000 random seeds),
1:10000,
# iteration code (formula ~)
~ {
set.seed(.x)
cysteine <- 0.3 * 10 ^ runif(1L, min = 0, max = 3)
adomed <- runif(1L, min = 0, max = 100)
setSpecies(
key = c("Cysteine", "adenosyl"),
initial.concentration = c(cysteine, adomed)
)
ss <- runSteadyState()
stopifnot(ss$result == "found")
list(
cysteine = cysteine,
adomed = adomed,
CGS = ss$reactions$concentration.flux[2],
TS = ss$reactions$concentration.flux[3]
)
}
)
# Combine all results and reshape the data
plotdata <-
scan %>%
bind_rows() %>%
gather("reaction", "flux", CGS, TS)
plot <-
ggplot(data = plotdata, aes(x = adomed, y = flux, group = reaction, tt_cys = cysteine)) +
geom_point(aes(color = reaction), alpha = 1 / 10, size = 3 / 4) +
labs(
x = paste0("Adomed (", getQuantityUnit(), ")"),
y = paste0("Flux (", getQuantityUnit(), " / ", getTimeUnit(), ")"),
color = "Species"
)
unloadModel()
plotly::ggplotly(plot, tooltip = c("tt_cys", "x", "y"))Parameter Estimation
This implements an example from the (Mendes2009 paper) on COPASI use cases.
loadSBML("http://www.ebi.ac.uk/biomodels-main/download?mid=BIOMD0000000010")
# get timeseries for the record
data_before <-
runTimeCourse(1000, 1)$result %>%
set_tidy_names(TRUE)
# read experimental data
data_experimental <-
read_tsv("data/MAPKdata.txt") %>%
rename(Time = time, "Mos-P" = "MAPKKK-P", "Erk2-P" = "MAPK-P") %>%
set_tidy_names(TRUE)
# define the experiments for copasi
fit_experiments <- defineExperiments(
data = data_experimental,
type = c("time", "dependent", "dependent"),
mapping = c(NA, "{[Mos-P]}", "{[Erk2-P]}"),
weight_method = "mean_square"
)
# define the parameters for copasi
fit_parameters <-
map(
parameter_strict(regex(c("V1$", "V2$", "V5$", "V6$", "V9$", "V10$"))),
~ {
val <- getParameters(.x)$value
defineParameter(parameter(.x, "Value"), lower.bound = val * 0.1, upper.bound = val * 1.9, start.value = val)
}
)
result <-
runParameterEstimation(
parameters = fit_parameters,
experiments = fit_experiments,
method = "LevenbergMarquardt",
updateModel = TRUE
)
# get timeseries for the record
data_after <-
runTimeCourse(1000, 1)$result %>%
set_tidy_names(TRUE)
plots <- list(
Erk2.P =
ggplot(mapping = aes(x = Time, y = Erk2.P)) +
geom_point(data = data_experimental, aes(color = "experimental")) +
geom_line(data = data_before, aes(color = "before")) +
geom_line(data = data_after, aes(color = "after")) +
scale_color_manual(values = c(before = "red", after = "green", experimental = "black")) +
labs(
x = paste0("Erk2-P (", getQuantityUnit(), ")"),
y = paste0("Time (", getTimeUnit(), ")"),
color = "Series"
),
Mos.P =
ggplot(mapping = aes(x = Time, y = Mos.P)) +
geom_point(data = data_experimental, aes(color = "experimental")) +
geom_line(data = data_before, aes(color = "before")) +
geom_line(data = data_after, aes(color = "after")) +
scale_color_manual(values = c(before = "red", after = "green", experimental = "black")) +
labs(
x = paste0("Mos-P (", getQuantityUnit(), ")"),
y = paste0("Time (", getTimeUnit(), ")"),
color = "Series"
)
)
unloadModel()
result$fitted.values
#> # A tibble: 2 x 5
#> fitted.value objective.value root.mean.square error.mean
#> <chr> <dbl> <dbl> <dbl>
#> 1 [Erk2-P] 10.08304 1.058460 0.3290446
#> 2 [Mos-P] 25.45444 1.681747 -0.2829374
#> # ... with 1 more variables: error.mean.std..deviation <dbl>
result$parameters
#> # A tibble: 6 x 8
#> parameter lower.bound start.value value
#> <chr> <dbl> <dbl> <dbl>
#> 1 (MAPKKK activation).V1 0.250 2.4643673 2.4643673
#> 2 (MAPKKK inactivation).V2 0.025 0.2465618 0.2465618
#> 3 (dephosphorylation of MAPKK-PP).V5 0.075 0.8817485 0.8817485
#> 4 (dephosphorylation of MAPKK-P).V6 0.075 1.4250000 1.4250000
#> 5 (dephosphorylation of MAPK-PP).V9 0.050 0.7280555 0.7280555
#> 6 (dephosphorylation of MAPK-P).V10 0.050 0.7070213 0.7070213
#> # ... with 4 more variables: upper.bound <dbl>, std..deviation <dbl>,
#> # coeff..of.variation.... <dbl>, gradient <dbl>
result$protocol
#> [1] "Algorithm started at 2017-09-06 22:19:03.\nFor more information about this method see: http://www.copasi.org/tiki-index.php?page=OD.Levenberg.Marquardt\n\nIteration 406: Objective function value and parameter change lower than tolerance (1/3). Resetting lambda.\n\nIteration 407: Objective function value and parameter change lower than tolerance (2/3). Resetting lambda.\n\nIteration 410: Objective function value and parameter change lower than tolerance (3/3). Terminating.\n\nAlgorithm reached the edge of the parameter domain 625 times.\n\nAlgorithm finished at 2017-09-06 22:19:17.\nTerminated after 411 of 2000 iterations.\n\n"
plotly::ggplotly(plots$Erk2.P)plotly::ggplotly(plots$Mos.P)